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X^ \ Abstract 

We describe a mechanism for pronounced biochemical oscillations, relevant to microscopic 
systems, such as the intracellular environment. This mechanism operates for reaction 
^ ■ schemes which, when modeled using deterministic rate equations, fail to exhibit oscillations 
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for any values of rate constants. The mechanism relies on amplification of the underlying 
stochasticity of reaction kinetics within a narrow window of frequencies. This amplification 
O ■ allows fluctuations to "beat the central limit theorem," having a dominant effect even though 

the number of molecules in the system is relatively large. The mechanism is quantitatively 
studied within simple models of self-regulatory gene expression, and glycolytic oscillations. 
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I. INTRODUCTION 

Biochemical oscillations have been studied from complementary experimental and the- 
oretical perspectives for many years. They appear to be generic processes in biological 
systems, with numerous examples known in both epigenetic and metabolic contexts [1-4]. 
Well known instances are circadian rhythms {e.g. in microorganisms [5,6]), and the oscilla- 
tion of ATP and ADP concentrations during phosphorylation of fructose-6-phosphate (F6P), 
a key step in glycolysis [7,8]. The study of biochemical oscillations has benefited from an 
intimate collaboration between experimentalists and theorists: reaction networks are pieced 
together in the laboratory, with hypotheses severely constrained by results from mathemat- 
ical models [4,8]. These mathematical models are generally constructed using deterministic 
chemical rate equations. Within this modeling framework, the signature of oscillatory be- 
havior is the existence of a limit cycle. The absence of a limit cycle is assumed to imply 
that the underlying reaction scheme, upon which the model is based, does not have cyclic 
behavior. Recently, a number of groups have challenged this assumption [9-12], mainly in 
the context of calcium oscillations. They have studied deterministic rate equations which 
support a limit cycle over a certain range of parameter values, and have found, from com- 
puter simulations, that the addition of noise can expand the region of parameter space in 
which cycles occur. These authors demonstrate that this effect can be maximized for par- 
ticular values of the system size or noise strength. On this basis, the effect has been termed 
"internal noise stochastic resonance." We show here that internal noise can have a far more 
profound effect on biochemical reaction kinetics; namely, inducing amplified oscillations in 
systems which, when modeled with rate equations, lack a limit cycle throughout their entire 
parameter space. Thus, from the view of conventional deterministic modeling, these reaction 
schemes would be immediately ruled out as candidates to describe biochemical oscillations. 

In this paper we develop a theoretical framework, using the Van Kampen system-size 
expansion [13-15], which allows an exact analytic description of these stochastically induced 
cycles - indeed, our theory reduces to a problem in linear algebra regardless of the reaction 
scheme under study, thereby allowing straightforward analysis. A crucial point is that 
this effect operates only in systems composed of a relatively modest number of molecules 
- typically in the range 10^ — 10^. Many biochemical reactions within cells operate with 
numbers of molecules in this range, which leads us to believe that intracellular processes 



may well exploit this amplification mechanism. Our result is at odds with the intuitive 
notion that fluctuations can be safely ignored for systems composed of many thousands of 
molecules - the reason being that, here, the amplitude of fluctuations is composed of two 
factors: the usual statistical factor 1/yN (where A^ is the typical number of molecules in 
the system) , and a large factor R^ 1 arising from the amplification of the underlying noise. 
If R/yN ~ 0(1) or larger, the intuitive conclusion that fiuctuations can be ignored for 
A^ ^ 1, is incorrect. 

For ease of presentation we shall introduce this mechanism through two relatively simple 
examples - one epigenetic (gene regulation), the other metabolic (glycolysis). Generaliza- 
tions to more complex reaction schemes is straightforward, since, as already mentioned, the 
analysis of fiuctuations reduces to an exercise in linear algebra. We will develop the the- 
oretical ideas using the language of the gene regulation model (defined in Figure la) for 
concreteness, because it involves only the populations of two constituents, unlike the model 
of glycolytic oscillations which involves four. 

The outline of the remainder of this paper is as follows. In section II we introduce the urn 
model representation of the gene regulation model, and proceed to instantiate this as a mas- 
ter equation. The deterministic limit of the master equation is derived in section III, and is 
shown to correspond to the usual chemical rate equation description. We calculate the effect 
of weak fiuctuations in section IV, using the Van Kampen system-size expansion. We obtain 
a description of the fiuctuations as a set of linear stochastic differential equations, which 
can therefore be analyzed exactly. The power-spectra obtained from these equations have 
resonance peaks which considerably enhance the 1/yN effects we would naively expect from 
fluctuations. The analysis (master equation, deterministic limit, flrst-order fluctuations) is 
repeated in a condensed form for Selkov's model of glycolysis in section V. Again, we flnd a 
strong peak in power spectrum of the fluctuations indicating amplifled oscillations. We end 
with a discussion of our results in the broader context of biochemical reaction networks. 

II. INDIVIDUAL-BASED STOCHASTIC MODEL FOR THE GENE REGULA- 
TION MODEL 

Consider, flrst, a simple reaction scheme in which an enzyme inhibits the transcription of 
its parent gene (Figures la and lb). More elaborate reaction schemes based on this simple 



model have long been proposed to explain circadian rhythms [5,16]. It is well-known that a 
deterministic model of this simplest reaction scheme involving two chemical agents (mRNA 
and the enzyme) will fail to produce cycles - at least three chemical agents (the mRNA, 
the enzyme, and an intermediate form, e.g. the primary RNA transcript) are required. 
It is important to recall that the deterministic model is an approximate representation of 
the actual chemical kinetics - it is strictly accurate for an infinite bath of molecules which 
are well-mixed. In order to describe a finite number of molecules it is necessary to use a 
discrete stochastic formulation in which one tracks the probability distribution of chemical 
concentrations over time. We have developed a standard stochastic treatment of the two- 
agent reaction scheme in Figures la and lb. The calculational steps used in the approach are 
illustrated in a flow-chart (Figure 2). By taking the limit of the number of molecules in our 
system, A^, to be infinite, we indeed recover the deterministic theory - this is an important 
benchmark. To handle the case of a finite number of molecules, we perform a system size 
expansion [13], which is a standard technique from the theory of stochastic processes, in 
which fluctuations are accounted for within a perturbative treatment. The second order 
terms in this expansion describe the fluctuations about the deterministic theory. These 
fluctuations satisfy linear equations and their statistics can be solved exactly. In particular, 
we calculate the power spectrum of these fluctuations. We find that for a wide range of 
parameters, the power spectrum has a pronounced maximum within a narrow window of 
frequencies. 

In order to perform systematic numerical comparisons between integrating the rate equa- 
tions and simulating the stochastic model it is convenient to envisage the reactions (Figures 
la and lb) occurring within two baths, as shown in figure 3a. The empty state has been 
replaced by null constituents Ei and E2 to be discussed further below. The reason for 
introducing two baths is simply because it is easier to set up the dynamics of the system 
and is closer in spirit to other systems (grounded in population dynamics) which have been 
analyzed in a similar way, and are well-understood [14,15]. 

Our aim is to build an individual-based stochastic model, and so the basic ingredients 
will be the number of M (mRNA) and Ei (null) constituents in bath 1 and of P (enzyme) 
and E2 (null) constituents in bath 2 — in a given stochastic realization at a given time 
t. The dynamics will consist of picking constituents from the baths at each time step and 
attempting the specified reactions. Performing many runs of such reactions will enable us to 



collect a large number of realizations and extract average behavior. The reactions proceed 
according to the rates shown in Figure 3a, and if the selection of constituents does not 
correspond to one of the four reactions, then the molecules are returned to their respective 
baths without any action being taken. 

On the basis of the specification so far, and summarized in Figure 3a, numerical simu- 
lations of the model can be carried out. We use the elegant Gillespie algorithm [17] which 
uses the information encoded in the reaction scheme to generate random time increments, 
in each of which a randomly selected reaction is forced to occur. The time increments and 
reactions are selected according to weighted distributions in such a way that the probability 
distribution of the stochastic time series generated is exact. The results of some of these 
simulations are displayed in Figures 4 and 5. 

It is also possible to derive a set of equations which describe the stochastic process which 
is defined by the model we have specified. To do this, we first note that there has been 
an implication that the choice of constituents at a given time step is a random process, 
only dependent on the numbers of the various molecules in the baths at that time, and not 
on choices or availability at previous time steps. Assumptions of this kind imply that the 
process is Markov, and so can be modeled using a master equation [13,18]. This is essentially 
a continuous time version of a Markov chain. Before we can write down this equation, we 
need to define some more quantities. 

Let us denote the number of molecules of the various kinds as follows: in bath 1 the 
number of molecules of M is rii and in bath 2 the number of molecules of P is n2. The state 
of the system is then denoted by the two numbers (ni,n2). We will frequently write this 
as n when we simply want to refer to the general state of the system. Note that we do not 
confer the status of independent variables on the numbers of Ei or E2 constituents. If we 
denote the total number of constituents in bath 1 by A^^i and that in bath 2 by A^2j then the 
number of Ei and E2 constituents are simply what is required to make up these numbers: 
{Ni — rii) El and (A^2 — 1^2) E2 constituents. This is why it was necessary to introduce the 
null constituents: if they had not been included, the number of M and P molecules would 
not have the freedom to vary — they provide room for independent changes in the numbers 
of both M and P molecules. 

We may now define transition rates from one state, n, to a different state n'. For instance, 
in the fourth reaction, the number of P molecules increases by 1 (recall that the number of E2 



molecules is not an independent variable). So in this case, n = {ni, n2) and n' = (ni, n2 + l). 
We denote the transition rate by T(n'|n). In our convention, initial states are on the right 
and final states on the left. When picking the constituents, there is a probability of rii/Ni 
of M being chosen, (A'"2 — n2)/A^2 that an E2 is chosen, and the reaction happens at a rate 
k2. This gives the result /c2(^i/^i)(^2 — 1^2) /N2 for this particular transition rate. Others 
may be found in the same way. A complete listing with n = {rii, ^2) is 



1. M^^Ei, n' = (ni-l,n2) 



2. P^^E2, n' = (ni,n2 



3. Ei^M, n' = (ni + l,n2). 



nri!\n) = lii^^. 



T(n1n) = /i^^ 



T{n \n) = fci = fcoexp (-An2/A^2) 77 



4. M + E2^M + P, n' = (ni,n2 + l). 



Tin \n) = k2- 



■Ni N2 



Note, we use an exponential function exp(— A[P]) to model the down-regulation of transcrip- 
tion. The primary reason is this form introduces one extra parameter (A) into the model, 
unlike the typical Hill form [16] (1 + q;[P]'')^^ which introduces two new parameters {a and 
q). Furthermore, in deterministic studies (albeit using more complex three-component mod- 
els), it is found that large values of q (~ 10) are required for limit cycles to occur [16]. For 
such large values of q, the Hill form is a rapidly decaying function, and it seems reasonable 
to replace it by a simple exponential function. 

Having defined the transition rates, we are now in a position to write down the master 
equation. It has the general form [13,18] 

j^P{n, t) = J2 T{n\r^)P{v!, t) - Y, T{n[\n)P{n, t) , (1) 

where P(n, t) is the probability that the system is in the state n at time t. This equation 
has a simple interpretation: the first term on the right hand side is the sum of the transition 
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rates into the state n from all other states n' and the second term is the sum of the transition 
rates out of the state n into all other states n'. When the second term is subtracted from 
the first, it gives the rate of change of the probability P. 

So far we have formulated a discrete stochastic model and given specified analytic forms 
for transition probabilities between states, which can be used in conjunction with the master 
equation ((T)) to, in principle, solve for the probability, P{n,t), that the system is in state n 
at time t. We will now start from the master equation and (i) determine the form that the 
model takes in the mean- field limit, and (ii) implement the system-size expansion to study 
the origin of the cycling behavior, which is absent in the deterministic (mean- field) limit. 

III. THE DETERMINISTIC LIMIT 

A straightforward way of obtaining the deterministic version of the model, valid in the 
limit of very large system size, from the master equation is to multiply (^ by rii and n2 in 
turn and to sum over all the states. This generates rate equations which are the deterministic 
equations if correlations between the variables are ignored. Let us illustrate the method in 
the case of ^2. We wish to calculate (^2) = J2n ''^^PiXLi t) by multiplying the master equation 
by n2 and summing over n. On the left-hand side we find d{n2)/dt. On the right-hand side 
are two terms which are nearly equal and opposite — if it was not for the n2 factor, they 
would be. The only difference between the two terms is that n and n/ are interchanged. So 
if n2 does not change in a reaction, as for reactions 1 and 3, then the two contributions do 
in fact cancel out. If it decreases by 1, as in reaction 2, then a shift in the sum over n2 gives 
an overall contribution of —1, and if it increases by 1, as in the reaction 4, then a shift in 
the sum over n2 the other way, gives an overall contribution of +1. This gives the result 

d , , / n2\ , / ni {N2 - n2) 



dV ' " \N2/ \Ni N2 

So far this is exact. The mean- field approximation enters through ignoring correlations, 
which vanish as Ni,N2 —>■ 00. So, for example, this means that (nirij) = {ni){nj) for 
i,j = 1, 2. If we make this approximation and introduce the fractions of M and P to be 0i 
and 02 respectively, in the limit Ni, N2 —^ cxd, then we find 

^ = -^02 + ^01(1-02). 
dt N2 N2^ ' 



This is the required deterministic equation. So in summary, if we define (pi = Ui/Ni and 
scale the time by introducing r = t/A^i, then the deterministic equations corresponding to 
the individual based stochastic model we have defined are 



d(t) 



dr 



- = -iii(j)i + koexp{-X(j)2)il-(j)i) , (2) 



O- ^ — = -/i202 + k2(pi (1 - 02) , (3) 

where a = N1/N2. Note that the mean-field approximation as applied to the first equation 
implies that (exp (— An2/A^2)) = exp (— A(r22)/A^2)- 

Most of the theoretical investigations of biochemical reactions start from a set of differ- 
ential equations of the type Q and (jS)). One of the first quantities of interest in such studies 
are the fixed points of the system. If we denote the fixed points with an asterisk and define 
X = A025 then X satisfies the transcendental equation 

7 -X _ Ail/^2-^ f.^ 



with the fixed points for (pi given by 



(5) 



^^ k2{X-X)' 

Note that 1 — 0]^ > and so the denominator of the right-hand side of Eq. (0} is never zero. 
Also, since the left-hand side of this equation is monotonically decreasing from its value at 
X = and the right-hand side is monotonically increasing from its value at X = 0, it follows 
that there is always a unique solution for X and therefore always just one fixed point. The 
condition < 02 — 1 implies that < X < A, but the condition < 0^ < 1 gives the 
stronger constraint 

X<-^^. (6) 

^2 + /i2 

The stability of these fixed points would always be of interest, but the question takes on 
an added significance in our case, since the stability matrix associated with the non-trivial 
fixed point (Eqs. (0)) and (0)) plays a central role in the analysis of the cycling phenomenon. 
To determine it, let us write the mean-field equations Q and Q in the form d(j)i/dT = fi{(j)) 
where i = 1,2. The fixed points are found from solving /j(0) = and a linear stability 
analysis consists of writing 0j = 0* + {(pi/ ^/Wi), where the (pi are small deviations from the 
fixed point [19] . Note that we have included extra factors of I/a/A^ in the small deviations 
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from the fixed point — which simply amounts to a re-scahng of 0j — compared to the usual 
linear stability analysis, so as to make contact with the system-size expansion in the next 
section. Linearizing about the fixed point gives d(f)i/dT = ^ . Mij(f)j, where Mij is defined 
by Mij = dfi/d(f)j\Fp. Here FP means "evaluated at the fixed point". The explicit forms for 
the entries of the matrix M are 

Mil = -/ii -/i:oexp(-A02) , 

Mi2 = -A;oA^'/'exp(-A0;)(l-0t) , 

M21 = k2a'/^ (1 - 0;) , 

M22 = -Crfi2 - Crk2(j)l . (7) 

These entries may be rewritten using the fixed-point equations to give 

Mil 



fll 


Mi2 = 
M22 = 


02 


1-01' 

,, ^1/2 *i^2 

01 



M21 = /i2al/2 5 , M22 = -^^2^ . (8) 

01 02 

From these expressions it is clear that all entries of M have a definite sign: Mii,Mi2 and 
M22 are negative, and M21 is positive. Therefore, the determinant and the trace of M are 
positive and negative respectively for all parameter choices. This already tells us that the 
fixed point is stable, but further analysis is required if we wish to know whether or not the 
the fixed point is approached in an oscillatory fashion. This will be discussed again in the 
next section. 

IV. ANALYSIS OF THE FLUCTUATIONS 

In contrast with the innate discreteness of the stochastic model, the mean-field equations 
involve functions of continuous variables. This is one of the reasons why they are more 
amenable to analytic treatment. It is the limit Ni,N2 -^ 00 which leads to the continuity 
of the mean-field equations, as well as to the elimination of the fluctuations. In the section 
we will discuss how we can keep continuity, but still not lose the stochastic nature of the 
system. This is achieved by using new continuous variables xi,a;2 in place of the previously 
used discrete variables ni,n2 to describe the probability distribution. The explicit form of 
the replacements are 
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The 1/yN terms are present since, by the central-hmit theorem, we expect fluctuations to be 
of the order of 1/yN when the variables ni, n2 are expressed in terms of the fractions ni/Ni 
and n2/N2- As Ni, N2 -^ 00, these fluctuations vanish, and the system is entirely described 
by the mean- field variables 0j(t), which can be found, in principle, by solving Eqs. (j21) and 
(jni) with given initial conditions. If we imagine a plot of the probability distribution P along 
the vertical axis, and the variables n^ along the horizontal axes, for various values of r, then 
at the initial time r = P is a delta- function spike at the starting values of the n^. As r 
increases, not only does the position of the peak move, but the probability distribution also 
broadens due to fluctuations. The (piij), which are the solutions of Eqs. (JD) and Q, tell us 
the position of the peak of the distribution and the variables a;j(r) tell us something about 
the distribution itself. Actually it turns out that to order 1/yN (here we use A^ to mean 
A^i or A'2, since they are assumed to be of the same order), the probability distribution 
is Gaussian, and since the position of the peak is specified, all that is left to determine 
is the variance. Higher terms in the 1/yN expansion give deviations from the Gaussian 
form, which our numerical simulations show are very small for reasonable values of A^. Once 
the leading (deterministic) contributions have been subtracted out from (jHl) (by, in effect, 
continuously moving the origin to the peak of the probability distribution) only fluctuations 
remain. If the l/v^Ai and 1/ \/N2 terms are then factored out, we may once again take 
N ^ 00. In this way the Xj become continuous variables, and terms of different orders can 
be identified in the master equation. 

The actual implementation of the system-size expansion is straightforward, if tedious. 
It is discussed clearly in [13] and in some detail in [14]. We shall illustrate it by applying 
it to one term in the master equation only. The reader should be able to understand the 
essential features of the method from this, and can then consult the above references to get 
a broader picture. Let us take reaction 2, in which the number of P molecules decreases 
by 1. It gives a contribution to the term T(n'|n)P(n, t) in the master equation (^ which 
is equal to /i2(^2/A^2)-P(^i5 n2,t). It also gives a contribution to T(?7.|n')P(n', t) in the same 
equation which is equal to /U2([^2 + 1]/A'2)-P(^i, ^2 + l,t). We can combine these two terms 
as follows: 



Reaction 2 : {S2 



f^2^fP{ni,n2,t) 

I\2 



(10) 



where S2 is a step operator which is defined in terms of its action on functions of the rii by 
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E^ "^'^{nx-.Ti'i) = ^{rii, 77-2 ± 1). Similar operators can be defined for tlie otlier variables. Tlie 
advantage of using these operators is that, within the replacement scheme 0, they have a 
simple form for large N. For example, 



S. 



±1 

2 



1± 



d 1 
+ 



d^ 



+ 



(11) 



^dx2 ' 2N2dxl 

Substituting Q and (fTTj) in (fTUI) . and expanding in inverse powers of y/N2 one finds a 
contribution proportional to l/v^A^: 



1 



/N, 



:fl2(f)2 



0X2 



and a term of order I/N2: 



N, 



1^2 






(12) 



(13) 



^2 ^ ^-^2. 

There are higher order terms, but this is as far as we need to go in the expansion. The 
quantity 11 which appears in Eqs. ()12p and (fT^ is numerically equal to P, but is instead a 
function of the Xj and of t. This means that the left-hand side of (^ now reads 
dP _ dH _ ^d(f)i dU _ /^rf02 dU 



dt dt dxi 



1 dli 

iVi 97 



Hsn 



a 



dt 8x2 

1 1 t/02 <9n 



(14) 



//Vj" (ir 9xi \/N2 dr 8x2 

Equating the right-hand side of the master equation (Eqs. (|T^ and (fT!?jl ) to the left-hand 
side (Eq. (HU) order by order, we find that a~^d(f)2/dT has a contribution —fi2(p2 (the 911/(9x2 
cancel out) and 



an 

a7 



Ai20- 






+ 



(15) 



where the dots mean that the reactions other than 2 will also give a contribution. 

This partial, but explicit, calculation allows us to see more clearly how the expansion 
works. At leading order we have a contribution to the mean-field equation for a^^d(f)2/dT 
which is equal to —fi24>2, which indeed appears in Eq. (jH)), and in none of the other mean-field 
equations. Therefore, an alternative, and in some sense more systematic, way of obtaining 
the mean-field equations is as the leading order terms in the large- A^ expansion. The next- 
to-leading terms give a partial differential equation for the probability distribution n(x, t) 
which, when we include all the other reactions, has the form 






d_ 
dx 






(16) 



«j 
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From p5|) we see that ^2(0;) = — /i2crx2 + • • . and Bu = yU2cr02 + . . .. In fact when all 

the reactions are included, the Ai{x) remain linear functions of the Xj and the Bij remain 

independent of them. We may therefore write 

2 
A{x) = J2Mi,x,. (17) 

i=i 

This means that the probability distribution at next-to-leading order, II{x,t), is completely 
determined by two 2x2 matrices: M and B, whose elements are independent of the Xj, 
and only functions of the (pj. In fact it is a characteristic of the large- iV expansion that M 
is nothing else but the Jacobian matrix with elements dfi/d(f)j, which may be calculated 
directly from the mean-field equations Q and Q- If the initial transients have died away 
and the solution to the deterministic equations has approached a fixed-point (p*, then M 
will simply be the stability matrix at that fixed point. The stability matrix for the fixed 
point Q and © is given in the previous section. The matrix B has entries 

Bn = /ii0i + /i;oexp(-A02)(l -0i) 

B22 = /U2(T02 + fc2Cr01 (1 - 02) 

Bu = B2i = 0. (18) 

So, to summarize the position so far, we know the Mij and B^j in terms of the parameters 
which define the individually-based stochastic model, and therefore by solving Eq. (fTT)|) we 
can find all we need to know about the fluctuations for large A^^i and N2. The partial 
differential equation (fTBj) is a Fokker-Planck equation — a continuous version of the master 
equation — and can be solved, in principle, given the initial condition that 11 is a delta- 
function spike at r = 0. In fact, for the simple case where the B^j are independent of Xj 
and the Ai are linear functions of the Xj, it can be solved exactly [13]. The result is a 
multi-variate Gaussian with (x,) = 0, as already mentioned. 

Our main aim in this paper is to understand oscillations and for this one of the main tools 
is Fourier analysis. The form of Eq. ()16p is not so useful for this purpose, but fortunately 
there is a completely equivalent formulation of the stochastic process which is ideally suited 
to investigation using Fourier transforms. Rather than write an equation for the probability 
distribution function 11, an equation for the actual stochastic variables Xi{T) can be given, 
in other words, the problem may be formulated as a set of stochastic differential equations 
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of the Langevin type. The Langevin equations which are equivalent to ()16|) are [13] 

^ = A{x) + r^,{T) , (19) 

where ?7j(t) is a Gaussian noise with zero mean and with a correlation function given by 

{v^irHir')) = B,,6iT - t') . (20) 

The system defined by Eqs. (fTn|l and pOjl is ideally suited to Fourier analysis, since the 
equations ^T^ are linear (since the Ai are) and Eq. ^F)^ implies that the noise is white, that 
is, the Fourier transform of its correlation function is frequency-independent. 
Taking the Fourier transform of (fT^ gives 

2 

-iuJXi{uj) = ^ MijXj{uj) + fji{uj) , (21) 

i=i 

where the tilde denotes the Fourier transform. We may write this as 

2 

y^ (^ij{uj)xj{ijj) = fji{ijj) , ^ij{oj) = -iuj5ij - Mij . (22) 

The Fourier transform of rii{t) has the correlation function 

{U^)fl,{uo')) = B,j{2'K)6{uo + uo') . (23) 

From Eq. (j^ we obtain Xi{uj) = ^ . ^^j^(Lu)fij{ij), and averaging the squared modulus of 
Xi gives the power-spectra 

2 2 

P.{u) = (|x.(a;)H =Y.Y. '^i/i^)B,k {<!>%! (uj) , (24) 

i=i fc=i 

where we have used ^ij{—iu) = ^j^{uj). [We have omitted the proportionality factor 2tt6{0). 
In practice, when comparing the analytically calculated power spectra to those generated 
from a numerical time series, one uses a discrete Fourier transform, and the proportionality 
factor is simply equal to the time increment of the recorded values in the time series.] 

To isolate the resonance, we note that Pi{uj) has the form of the ratio of two power series 
in u). The denominator, which will largely control the position of the resonance, can be 
simply expressed as the determinant of the matrix $. If we define D{uj) = det^{uj), then 
the denominator is just \D{uj)\'^. 
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In Eqs. fj21|) - ()24|) we have written everything in a rather general form, since the forniahsm 
can be seen to generahze to the case of an arbitrary number of constituents rather easily, 
such as in the model of glycolysis in the next section which has 4 constituents. However, for 
the model we are presently considering, there are only two constituents and the expressions 
for Pi{uj), i = 1,2 have simple explicit forms: 



\D(uj) 



P^ico) = mco)\') = ^^^^r^ , (25) 



where D{lu) = —ur + iuj trM + det M and where 



02 = BxxMl^ + B22Ml^ , /?2 = ^22 . (26) 

Since the elements of the M and B matrices are known in terms of the parameters of the 
model and the fixed point values of the 0j, so are the ctj and /3j. 
The denominator of the power spectrum in Eq. (^3]) is given by 

2 , /4- A#\2 2 



\D{yj)f = (cj2-detM) +(trM) 



(jj 



,2 



= {iu'-niy + T'uj\ (27) 

where, since in Section lTTTl we found det M > and trM < 0, we have introduced Qq = det M 
and r = —trM. So, we may write the power spectrum as 

R^ico) = (|x.(c.)H = ^ ^^±^ ^ . (28) 



UJ 



2 - fi2)2 + Y^U^ 



This form for the power-spectrum shows clearly the existence of a resonance: for a value of 
cu^ the denominator becomes small, and the power spectrum has a large peak centered on 
this frequency. 

To analyze the nature of the resonance in more detail, let us set z = uj'^ and ask for what 
values of z the power spectrum (j^ has a maximum. The condition dP/dz = gives 

(3z^ + 2az + [a {T^ - 2^^) _ /j^^j _ q ^ ^29) 

where we have dropped the index i on Oj and Pi. Let us begin by neglecting the term Piu'^ 
in the numerator of Eq. ()28|) which may be justified for some parameter choices. Then the 
condition (|^ simply becomes z = {2Ql — T^)/2. Since we require z = cu'^ > 0, this implies 
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2Ql > r^. In terms of the stability matrix M, this condition reads 2detM > (trM)^, 
which imphes that the eigenvalues of M are complex. In the last section we showed that the 
eigenvalues of the stability matrix were either real and negative or complex with a negative 
real part. The condition that the power spectrum has an extremum imposes the latter. 

In reality, the presence of the factor Piu"^ in the numerator of the power spectrum may 
have a significant effect, and the full condition (|^ has to be used. From Eqs. (jTHj) and (|^ 
we see that the a, and jSi are always positive, and so from Eq. (j^^ we see that the sum of 
the roots of this equation is negative and so at least one of the roots is negative. For the 
other one to be positive we require 



a 



(r^ - 2nl) -/3n^Q<o. (30) 



This is the general condition for an extremum of the power spectrum to occur. It goes 
beyond the the previous condition, which only involved the eigenvalues of the matrix M, 
since it contains the (3i which come from consideration of the fluctuations about the fixed 
point. We may now ask if the extremum is a maximum. It is straightforward to show by 
calculating dPP/dz^ that, if a positive solution to Eq. ()29|) exists, then it is automatically a 
maximum. 

A peak in the power spectrum indicates a resonant frequency, and a corresponding oscil- 
lation at this frequency. [We must stress that our use of the word "resonance" here follows 
the more general usage; namely, a peak in the power spectrum as a function of frequency. 
In "stochastic resonance" and its derivatives [9,11] (and references therein), the "resonance" 
refers to a maximized response of the system as a function of some parameter, such as the 
strength of the internal noise. Because of the potential for confusion, we often use the term 
"amplification" for the effect under consideration in this work.] The height of the peak 
indicates the strength of the amplification. The width of the peak indicates the amount of 
frequency dispersion one would observe about the resonant frequency. We have calculated 
the power spectra exactly for this epigenetic model. For a specific set of parameters, we show 
the power spectra for fluctuations in mRNA and enzyme concentrations (Figures 4c and 4d 
respectively) . We have also computed stochastic realizations of time series for this reaction, 
using the exact Gillespie algorithm [17] (Figures 4a, b). Note, the resonant amplification 
R for mRNA (defined as the height of the peak of the power spectrum relative to -P(O)) is 
quite pronounced {R ~ 30), and, indeed, regular oscillatory peaks are easily discernible in 
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the corresponding time series (Figure 4a). The constant dashed hnes in Figures 4a, b are 
the results from the deterministic theory - which predicts a complete absence of cycling. We 
emphasize that there are two competing quantities which combine to determine the impor- 
tance of these amplified oscillations - these are i) the relative height of the peak in the power 
spectrum R, and ii) the mean number of molecules in the system A^. The amplitude of the 
oscillations in a time series measurement of the concentration will be of the order R/yN . 
This is why the effect dies away for macroscopic systems in which A^ is extremely large. In 
the example given here i? ~ 30 (for mRNA fluctuations), and so we can estimate that the 
oscillations will be negligible in intracellular systems for which A'^mRNA ^ 1000. Given the 
simplicity of this two component model, and the modest estimate above, this mechanism 
may be relevant to circadian rhythms in procaryotic cells, such as cyanobacteria [5], since 
such cells are small, and have no nuclear membrane separating transcription and translation 
processes (justifying, to some degree, a two-component model). 

In Figure 5 we compare the exact form for the mRNA power spectrum, given above 
in ()25|1 . to a numerically generated power spectrum. The latter is computed from discrete 
Fourier transforms of a large number (10,000) of long (60,000 iterations) time series generated 
using the exact Gillespie algorithm described earlier. There is excellent agreement between 
the exact and numerical spectra, as expected. 

V. SELKOV'S MODEL OF GLYCOLYSIS 

We now turn to an example from intracellular metabolic dynamics. Consider the key 
step in glycolysis in which the enzyme PFKl catalyzes the phosphorylation of F6P. This 
reaction is actually rather complex, and rather sophisticated models of this process have been 
developed over the years, guided by ever more quantitatively precise experiments [3,7,8]. 
For the purposes of illustration, we shall concentrate on one of the first models of this 
reaction due to Sel'kov [20] . Although this model is no longer regarded as being an accurate 
representation of this reaction, it captures the essence of the process. Sel'kov's reaction 
scheme is illustrated in Figures Ic and Id. A key parameter in this scheme is the number 
of ADP molecules (conventionally denoted by 7) required to activate the PFKl enzyme. 
Within the deterministic modeling framework it has been shown that 7 > 1 is required for 
cyclic behavior to emerge from the model (even though the biochemistry of PFKl demands 
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7 = 1), and that even then, the cychng exists in a very narrow range of parameter values. 
We have reformulated Sel'kov's reaction scheme using the stochastic framework, but have 
restricted our attention to the more biologically plausible case of 7 = 1. 

As before, the deterministic model is retrieved when the number of molecules in the 
system is taken to be infinitely large, and predicts constant, non-cycling, concentrations. 
When this number is finite, we account for the fluctuations in the system using the system 
size expansion, and solve the simple linear theory in order to calculate the power spectra for 
the various chemical agents. 

The reactions which comprise the Selkov model of a key stage (phosphorylation of F6P 
by the PFKl enzyme) in glycolysis are shown schematically in Figures Ic and Id. Writing 
them out in a slightly different format, they read 

ADP + PFKl ^ PFKl/ADP 

k-3 

^^ATP 
ATP + PFKl/ADP ^ PFKl/ADP/ATP -^ PFKl/ADP + ADP 

k-i 

ADP^^ 

Once again we will introduce two baths. 

For notational simplicity we will denote ATP by 5*1, ADP by 5*2, PFKl by E2, PFKl/ADP 
by A and PFKl/ADP/ATP by B. The first two will be placed in bath number 1, and the 
last three in bath number 2 (Figure 3b). It will also be necessary to introduce "nulls" into 
bath 1, which we will denote by Ei. These represent the aqueous space which can potentially 
be filled by ATP or ADP. With this notation A is defined to be E2S2 and B is defined to 
be S1E2S2 and the five reactions are 

1. E, ^^ S, 



2. 


S2^E, 


3. 


Si + A^ B + Ei 

k-i 


4. 


B + Ei^A + S2 


5. 


S2 + E2 ^ A + Ei 

k--i 
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Let us denote the number of molecules of the various kinds as follows: in bath 1, the number 
of molecules of 6*0,, a = 1, 2 is nia, and in bath 2, the number of molecules of A is rii and of 
B is n2. The state of the system is then denoted by the four numbers {mi,m2, ni, n2), which 
we will again write as n when we simply want to refer to the general state of the system. If, 
as before, we denote the total number of constituents in bath 1 by A'^i and that in bath 2 by 
N2, then the number of Ei and E2 constituents is simply Ni — rrii — 1712 and N2 = rii — n2, 
respectively. 

The transition rates for this model are: 

1. n' = (mi + 1, m2, ni, 722). 

(A^i - mi - m2) 



T{nf\n) = Ui- 



2. n' = {mi, 7712 - '^,ni,n2). 

Tin \n) = ^^2^ • 

3. Forward reaction: n' = (mi — 1, ?Ti2, rii — 1, 722 + 1). 

^(^"^) = ^^i^iV2- 

Backward reaction: n' = (mi + 1, m2, ni + 1, n2 — 1). 

^(^1^) = ^-^ Ni iV2 

4. n' = (mi, m2 + 1, ni + 1, n2 - 1). 

(A^i - mi - m2) n2 



Tfn'ln) = k 



2" 



A^i Ar2 

5. Forward reaction: n' = (mi, 7712 — 1, ni + 1, 712). 

7712 {N2 -7li- 7I2) 



T(n'|n) = A;3 



A^i N2 



Backward reaction: n' = (mi, m2 + 1, ni — 1, n2). 

rr( n \ 7 (^1 - "^1 - "^2) ni 
T(n n = A;_s 
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To find the deterministic equation we define Sa = rria/Ni, a = ni/N2 and b = n2/A^2 and 
scale the time by introducing r = t/Ni . This gives the deterministic equations, correspond- 
ing to the individual based stochastic model we have defined, to be 

— — = z/i(l - si - S2) - hsia + A;„i(l - si - S2)b , (31) 

ar 

ds 

—— = -V2S2 + ^2(1 - si- S2)b - ksS2{l - a-b) + A;-3(l - si - ^2)0 , (32) 

ar 

a~^ -— = -kiSia + A;_i(l - Si - 52)^ + ^2(1 - Si - 52)6 
ar 

+^3^2(1 - a-b) - A;-3(l - si - S2)a , (33) 

(7~^—- = kiSia - k_i{l - si - S2)b - k2{l - Si - 52)6 , (34) 

UT 

where a = Ni/N2. The equations (j!?T|l - (|H^ are Selkov's equations [8,20], except for the 
additional factors of (1 — si — S2) which arise from the introduction of the nulls Ei. 

To investigate the fixed point structure, let us first note that if we replace the term 
(1 — si — S2) by Q in the equations (j^ - ljH^ . it will enable us to examine the Selkov model 
{Q = 1) and our deterministic model fi = 1 — si — S2, in tandem. From Eqs. (j31|) and (|34p 
we see that the fixed point has either to have fi* = or z/i = k2b* (all fixed point values are 
again denoted by asterisks). The first condition cannot hold in the original Selkov model, 
but can in our version of the model: in this case we see from Eq. (jH^ that Sg = 0, and so 
s* = 1. This implies a* = and b* is indeterminate. If f2* 7^ 0, remarkably a unique fixed 
point is found in both forms of the model, and moreover it takes on a reasonable simple 
form. The specific results are: 

• Fixed points of the original Selkov model 

W2k-3 + U^k3] f z/iV^ , l^i 

(fc2 + fc-i) 1-1- , S2 = — , 



kik2ks V k2J 1^2 



/■'=3 (i-^V 6- = ^. (35) 

[i'2k-3 + uik^] \ k2j fc2 



Fixed points of the modified Selkov model 

The trivial fixed point: s\ = l, S2 = , a* = 1. 
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The non-trivial fixed point: 

W2k-3 + i^ih] {k2 + k^i] 



A 



So = — 



[z/2A;_3 + uiks] V k 



i'2 



J^l 


kik2k3 


( 

1- 

V 


y{\ 


1^2 


A 


k2j 


J^l 








k2 


•) 







b* = j^, (36) 



where 



A = kik2k3 (l - ^) (l + 7) + h/c-s + Jyiks] {k2 + k-i) . 

To make contact with the discussion in earher sections, let us introduce the notation 
01 = Si,02 = -52,03 = a, 04 = b. Tlicu the mean-field equations (ISII)-(jSl|) take the form 
d(f)i/dT = fi{(j)) where i = 1, ... ,4. Performing a linear stability analysis gives the matrix 
Mij. The explicit forms for the entries of this matrix for the case when 0* is the fixed point 
(J36|) are given in the Appendix. 

To study the fiuctuations we introduce new continuous variables xi,X2,X3, x^ in place of 
the previously used discrete variables mi,m2,ni,n2, to describe the probability distribution. 
The explicit form of the replacements are 

mi , , 2;i "^2 , , 3^2 

n + -7^^ 5 T^ = 02 + 



A^i ViVi A^i "^ v^' 

— = 03 + -^ , — = 04 + — ^ . (37) 



Proceeding as in section HVl we carry out a system-size expansion. We find the deterministic 
equations (jHT|) -()34 |) to leading order and the Fokker-Planck equation ()16p at next-to- leading 
order, with Ai{x) given by Eq. |T7jl . except now that j runs from 1 to 4. So in this case, 
the probability distribution at next-to- leading order, n(a;, r), is completely determined by 
two 4x4 matrices M and B. The stability matrix for the fixed point ()3fj|l is given in the 
Appendix, where we also give the entries of the matrix B. The analysis of the Fokker-Planck 
equation, the conversion to a set of Langevin equations and the subsequent Fourier analysis 
is now exactly as in Section IIVI but with z, j = 1, . . . , 4. 

We have calculated the power spectra for the ATP and ADP concentration fiuctuations 
exactly. We find that these power spectra have very large and sharp peaks for a wide range of 
parameter values (Figures 6c and 6d) - indicating that the concentrations of these molecules 
will undergo violent oscillatory behavior within a small-system setting. Explicit stochastic 
simulations of the reaction network show that indeed large amplitude cycling occurs (Figures 
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6a and 6b). The constant dashed hnes are the predictions from the corresponding deter- 
ministic theory - showing a complete absence of cychng, as expected. The power spectrum 
peaks found in this example are rather large - of the order of 150 for ADP. This indicates 
that intracellular environments composed of up to 25,000 ADP molecules operating accord- 
ing to this reaction scheme will show pronounced oscillations, caused by amplification of the 
underlying stochasticity of the reaction kinetics. 

VI. DISCUSSION AND CONCLUSIONS 

We have described a mechanism by which a microscopic biochemical system can, loosely 
speaking, "beat the central limit theorem" via an amplification of stochastic fiuctuations in 
the reaction kinetics of its constituents. This mechanism leads to sizable oscillations in the 
concentrations of the reagents for reaction schemes, which when modeled using rate equa- 
tions, show no cycling behavior for any values of their rate constants. We have illustrated 
the effect in two very different biological examples - self-regulation of a gene, relevant to 
the study of circadian rhythms, and the dynamics of ADP, ATP, and PFKl concentrations 
during glycolysis. 

We stress here that the system-size expansion always leads to linear equations for the 
fluctuations, with coefficients related to the steady-state concentrations predicted from the 
flrst-order theory (i.e. the deterministic rate equations). Thus, the evaluation of the power 
spectra is simply an exercise in linear algebra. The signature of the ampliflcation mechanism 
is the existence of a peak in the power spectra. The existence of a peak is guaranteed if the 
stability matrix of the deterministic rate equations has complex eigenvalues, with negative 
real parts to ensure stability. This provides a simple test for the existence of the mechanism 
directly from the deterministic theory. Crudely speaking, if the approach to the deterministic 
steady-state occurs via damped oscillations, then the inclusion of second-order fluctuations 
will lead to the ampliflcation of sustained oscillations. It is also important to point out that 
the mechanism described here requires no external tuning of rate constants. This is because 
the underlying stochasticity has a flat spectrum in frequency space (i.e. white noise), and 
this automatically excites the resonant frequencies of the system. The oscillations are also 
robust - in the two examples given in this paper, cycles are present over a very broad range 
of parameter values. 
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The examples of circadian rhythms and glycolysis given here are very simple and designed 
to illustrate the amplification mechanism in two well-known areas of cell biology. It will be 
important to analyze more realistic (and hence, more complex) reaction networks using the 
same techniques. There is intense current interest in analyzing the dynamical properties of 
such networks, especially with regard to robustness [21]. Deterministic rate equations are 
used for these analyzes, but this may need to be revisited on the basis of the results of this 
paper. For a network involving M different chemical reagents, the analysis of fiuctuations in 
the system size expansion reduces to linear algebra with matrices of rank M . If oscillatory 
behavior is found, the "special frequencies" will be closely related to the M normal modes 
of the system. This implies an interesting analogy with the theory of small oscillations in 
mechanics [22] , whereby complex oscillatory motion can be decomposed into normal modes, 
with the lowest frequency modes almost completely characterizing the dynamics. As we 
have stressed, this mechanism only operates within a microscopic system, such as a cell, or 
internal region of a cell. There is, however, a means by which this mechanism can lead to 
macroscopic oscillations in a population of cells; namely, synchronization of the individual 
cellular cycles through weak intercellular interactions. 

The amplification mechanism described here should be considered as a possible cause 
for oscillatory phenomena in more complicated biological situations. The hypothesis is 
theoretically elegant in that it does not require ad hoc nonlinearities postulated simply to 
"generate" cycles. For example, endogenous circadian rhythms in eucaryotes appear to be 
controlled by a complex transcriptional feedback loop involving the genes dock, cry, Bmall 
and members of the per gene family [23,24]. Oscillations in the activities of these genes may 
be explained by this new mechanism if the number of mRNAs, translational complexes, and 
product proteins involved are small enough to allow the amplification term R to dominate. 
This is not unreasonable given that the number of mRNAs (for a given gene) in a eucaryotic 
cell ranges from 10^ — 10"^ [25]. This amplification mechanism may also explain oscillations 
within the more complex reality of glycolysis. Intracellular F6P concentration in the yeast 
Saccharomyces cerevisiae is approximately 0.1 mM [26]. Therefore, a cell of 10 /im diameter 
contains on the order of 10^ to 10^ F6P molecules. In mammalian cells, the concentration 
of PFKl appears to vary between 0.1 and 1 micromolar [27], yielding on the order of 10^ to 
10^ PFK tetramers per cell. Most eucaryotic cells have on the order of 10^ ATP molecules, 
but some orders of magnitude less ATP participates directly in glycolysis. Therefore, the 
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system size of glycolysis, especially if compartmentalized in a eucaryotic cell, approaches 
that required for noise amplification to have a pronounced effect. 

The condition for this effect of amplified oscillations is simply the existence of damped 
oscillations in the corresponding rate equation model {i.e. complex eigenvalues with nega- 
tive real part in the stability matrix). Thus, this mechanism is likely to be widespread in 
biochemical networks, implying that dynamics in the cytoplasmic environment may be far 
richer than formerly supposed on the basis of chemical rate equations. 



We thank Stuart Lindsay for introducing us to the subject of circadian rhythms in procary- 
otes. We acknowledge the NSF and NIH for partial support, under grants DEB-0328267 
and DMS/NIGMS-0342388. 
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APPENDIX A: EXPLICIT FORMS OF THE MATRICES M AND B 

As explained in Section of the main text, the matrix M, for the set of deterministic 
equations dcpi/dr = fi^cp), is defined by Mij = dfi/dipjlpp. For the system (| ^ -(|3 ^ the 
entries of the matrix are 

Mil = —^1 — kia — k^ib , M12 = —Vi — k^ib , 
Mi3 = -a^/^ [kisi] , Mi4 = a^/=^ [nk.i] , 

M21 = -k2b - k^sa , M22 = -{1^2 + ks) + (^3 - k^3)a + (k^ - ^2)^ , 
M23 = ai/2 ^^_^^ + ^^^^] ^ jvf24 = ori/2 [A;2fi + A;3S2] , 

M31 = (7^/2 p_^ _ ^^)^ _ (^^ + ^_^)^] ^ 

M32 = (T^/2 [k3 + {k_3 - k3)a - (A;_i + /C2 + fc3)fe] , 

M33 = a [-kisi - ^352 - k^^Q] , M34 = a [(A;_i + ^2)!! - ^352] , 

M41 = a^/2 [^^^ ^ ^_^^ ^ j^^^ ^ M^^ ^ ^1/2 p_^ ^ ^^)^| ^ 

M43 = a [kiSi] , M44 = CT [-{k^i + k2)n] , 

(Al) 

where Q = 1 — si — S2- The a, b, Si and S2 in these matrix elements are all assumed to be 
evaluated at the fixed point (jHSJ- 
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The noise-correlation matrix, Bij, is symmetric and given by 

Bii = z/i(l - si - S2) + /cisifl + A;_i(l - si - 82)6, 

B12 = 0, 5i3 = a^/^ [kisia + A;_i(l - Si - Sa)^] , Bu = -B13 , 

522 = l^2S2 + ^2(1 - Si - S2)b + k3S2{l - a - b) + k^s{l - Si - S2)a , 

B23 = o-^/^ [k2{l - Si - S2)b - k3S2{l - a-b) - k^^il - Si - 52)0] , 
^24 = a^/' [-A;2(l - si - S2)6] , 

-B33 = 0- [/ciSia + /i;_i(l - si - 52)^ + ^2(1 - si - S2)&] • 

+A;3S2(1 - a-b) + k-3{l - si - 52)0] , 
-B34 = 0" [-fciSia - A;_i(l - si - 82)6 - k2{l - si - S2)fe] , 

^44 = a [/cisia + A;_i(l - si - 52)^ + ^^2(1 - si - 52)^] • (A2) 
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Figure Captions 

Figure 1: Reaction networks for tlie two models considered in tliis paper. In tlie gene reg- 
ulation model (panels a,b), mRNA (M) and the corresponding enzyme (P) have rates of 
degradation /ii and fi2 respectively. Transcription occurs with rate ki (with inhibition), 
and translation occurs with rate k2- The functional form of inhibition is chosen to be an 
exponential function of relative enzyme concentration, with "regulation parameter" A. In 
Sel'kov's model for phosphorylation of fructose-6-phosphate (panels c,d), ATP enters the 
system with rate ui and ADP leaves the system with rate 1^2 ■ The enzyme PFKl is acti- 
vated by ADP, and this complex catalyzes the reaction with F6P (not shown) and ATP. 
Crucially, we assume that only one molecule of ADP is required to activate PFKl. 

Figure 2: A flow diagram indicating the connections between deterministic and stochastic 

models of the same reaction scheme. Within the deterministic model, if limit cycles are 

absent, then the model predicts that cycling is absent. The stochastic model is analyzed 

using the system size expansion. At leading order the deterministic theory is retrieved. The 

corrections account for stochastic fluctuations. If these are amplified (with a factor R), 

then the model predicts large cycles, so long as the typical number of molecules A^ satisfies 
ArV2 _ j^ 

Figure 3: In panel (a), the gene regulation reaction network (Figures la and lb) is recast 
in terms of two baths. Available units of space (or resources) for mRNA in bath 1 are 
represented by constituents Ei. Likewise, available units for the enzyme in bath 2 are 
represented by constituents E2. Selkov's model for phosphorylation of F6P (Figures Ic and 
Id) is similarly recast in terms of two baths in panel (b). ATP (Si) and ADP (6*2) molecules, 
along with available units Ei are contained in bath 1, while the PFKl enzyme {E2), and its 
complexes {A and B) are contained in bath 2. 

Figure 4: Dynamics of mRNA and enzyme concentrations in the simple gene regulation 
model defined in Figs, la and lb. The rate constants are ko = 1.0, k2 = 0.1, ^i = fi2 = 0.001, 
the regulation parameter is A = 100.0, and the bath sizes are Ni = N2/2 = 20480. Panels 
(a) and (b) show time series of the relative concentrations of mRNA {m) and enzyme {p) 
respectively. The constant dashed lines are the predictions from the deterministic theory. 
Panels (c) and (d) show the power spectra (normalized so that P(0) = 1) associated with 
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the fluctuations in the mRNA and enzyme concentrations respectively. The resonance peaks 
are the signature of this cycling mechanism. The relative amplification of fluctuations is 
R = 27.9 for mRNA, and R = 3.4 for the enzyme. 

Figure 5: The unnormalized power spectrum for mRNA fluctuations, as a function of fre- 
quency (for parameter values, refer to Figure 4c). The solid line is the exact result ()25|). while 
the circles correspond to the power spectrum computed from the discrete Fourier transforms 
of 10,000 mRNA time series generated from Gillespie's algorithm. 

Figure 6: Dynamics of ATP and ADP concentrations in the stochastic reformulation of 
Sel'kov's reaction scheme. The rate constants are ki = 1.0, k_i = 0.5, A;2 = 0.2, A;3 = 
0.2, k-3 = 1.0, z/i = 0.0005,1/2 = 0.03, and the bath sizes are A^i = A^2 = 4096. Panels 
(a) and (b) show time series of the relative concentrations of ATP (si) and ADP (S2) 
respectively. The constant dashed lines are the predictions from the deterministic theory. 
Panels (c) and (d) show the normalized power spectra associated with the fluctuations in 
the ATP and ADP concentrations respectively. The sharp resonance peaks are the signature 
of this cycling mechanism. The relative amplification of fluctuations is R = 15.1 for ATP, 
and R = 150.4 for ADP. 
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